#include "Solution.h"
using namespace blitz;

bool Solution::Analysis()
{
	Initial();
	return true;
}

bool Solution::Initial()
{
	const int edge_num = Object->GetEdgeNum();
	Z.resize(emx::Itmax);
	I.resize(emx::Itmax);
	V.resize(emx::Itmax);
	for (int i=0; i<emx::Itmax; i++)
	{
		Z(i).resize(edge_num, edge_num);
		I(i).resize(edge_num);
		V(i).resize(edge_num);
	}
	return true;
}

emx::real_data Solution::Singular(int patchIdx)
{
	emx::Real3D r[3], rho[3];
	emx::real_data p[3], sinBeta[3], _sinBeta[3], cosBeta[3], _cosBeta[3];
	emx::real_data result = 0.0;
	for (int i=0; i<3; i++)
	{
		r[i] = Object->GetPatch(patchIdx)->Vertex((i+1)%3)->Position - Object->GetPatch(patchIdx)->Vertex(i)->Position;
		rho[i] = Object->GetPatch(patchIdx)->Vertex(i)->Position - Object->GetPatch(patchIdx)->Center;
		sinBeta[i] = -dot(rho[i],r[i])/(sqrt(dot(rho[i],rho[i])*dot(r[i],r[i])));
		cosBeta[i] = sqrt(1-sinBeta[i]*sinBeta[i]);
		p[i] = cosBeta[i]*sqrt(dot(rho[i],rho[i]));
		_sinBeta[i] = dot(rho[(i+1)%3],r[i])/(sqrt(dot(rho[(i+1)%3],rho[(i+1)%3])*dot(r[i],r[i])));
		_cosBeta[i] = sqrt(1-_sinBeta[i]*_sinBeta[i]);
	}
	for (int i=0; i<3; i++)
	{
		result += p[i]*log(cosBeta[i]*(1+_sinBeta[i])/(_cosBeta[i]*(1-sinBeta[i])));
	}

	return result;
}